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BOUNDARY  CONDITIONS  IN  THE 
NAVY  COASTAL  OCEAN  MODEL 


1.  INTRODUCTION 

The  Navy  Coastal  Ocean  Model  (NCOM)  (Martin  2000)  is  an  ocean  general  circulation  model 
(OGCM)  being  developed  as  the  ocean  component  of  the  Coupled  Ocean- Atmosphere  Mesoscale 
Prediction  System  (COAMPS).  CO  AMPS  is  being  developed  for  operational  use  by  the  Navy 
(Hodur  1997).  In  addition  to  being  used  as  an  ocean  component  of  COAMPS,  the  NCOM  can 
be  used  as  a  stand-alone  OGCM  for  studies  of  ocean  dynamics  ranging  from  the  littoral  scale  of 
coastal  regions  up  to  basin  and  global  scales. 

NCOM  has  a  free  surface  and  is  based  on  the  primitive  equations  and  the  hydrostatic,  Boussi- 
nesq,  and  incompressible  approximations.  The  model  uses  an  Arakawa  C  grid  and  is  leapfrog  in 
time  with  an  Asselin  filter  to  suppress  timesplitting.  The  propagation  of  surface  waves  and  vertical 
diffusion  are  treated  implicitly.  A  choice  of  the  Mellor-Yamada  Level  2  or  Level  2.5  turbulence  mod¬ 
els  is  provided  for  the  parameterization  of  vertical  mixing.  The  horizontal  grid  is  curvilinear.  The 
vertical  grid  uses  sigma  coordinates  for  the  upper  layers  and  j2:-level  (constant  depth)  coordinates 
for  the  lower  layers,  and  the  depth  at  which  the  model  changes  from  sigma  to  z-\evel  coordinates 
can  be  specified  by  the  user.  A  source  term  in  the  model  equations  can  be  used  to  specify  input  of 
river  and  runoff  inflows. 

Boundary  conditions  (BCs)  have  been  implemented  in  the  NCOM  to  provide  forcing  at  the 
surface  and  at  lateral  open  boundaries.  At  the  surface,  the  BCs  that  can  be  applied  include 
surface  atmospheric  pressure  and  wind  stress,  penetrating  solar  radiation,  and  fluxes  of  scalar  fields 
including  heat  and  moisture  fluxes  for  the  temperature  and  salinity.  Surface  values  of  scalar  fields 
can  also  be  relaxed  towards  prescribed  values.  Local  tidal  forcing  can  be  applied  by  specifying  the 
horizontal  gradient  of  the  surface  elevation  due  to  the  local  tidal  potential. 

At  lateral  open  boundaries,  a  variety  of  BCs  can  be  used  to  specify  the  flow  into  and  out  of 
the  model  domain,  to  specify  the  values  of  scalar  fields,  to  provide  tidal  forcing,  and  to  allow  for 
gravity  wave  transients  and  disturbances  to  propagate  out  of  the  region.  Most  of  these  BCs  have 
been  implemented  as  “tendencies”  towards  desired  boundary  values  with  time  scales  appropriate 
to  the  particular  phenomena. 

The  NCOM  also  provides  for  the  use  of  cyclic  lateral  BCs  in  either  or  both  horizontal  di¬ 
rections.  Besides  allowing  for  periodic  domains,  these  BCs  also  allow  the  NCOM  to  be  run  as  a 
two-dimensional  (2-D)  vertical  section  or  as  a  one-dimensional  (1-D)  profile  with  variability  only 
in  the  vertical,  i.e. ,  like  a  local  mixed-layer  model.  The  2-D  and  1-D  capabilities  are  very  useful 
for  running  test  problems. 

Manuscript  approved  January  30,  2001. 
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In  this  report  we  document  the  BCs  that  have  been  implemented  in  NCOM  version  1.5.  Sec¬ 
tion  2  introduces  some  notation  that  will  be  used  throughout  the  remainder  of  the  report.  Section  3 
describes  the  cyclic  BCs  available  and  how  to  configure  the  NCOM  for  1-D  or  2-D  simulations. 
Section  4  presents  the  surface  forcing  that  can  be  applied  at  the  air-sea  interface  with  particular 
attention  to  the  required  physical  units.  Section  5  discusses  the  surface  BC  for  the  turbulence 
submodels.  Section  6  explains  the  various  possible  BCs  that  can  be  applied  along  the  NCOM  open 
boundaries  for  the  elevation,  barotropic  and  baroclinic  velocities,  vertical  velocities,  scalar  fields, 
and  the  vertical  mixing  variables.  Section  7  discusses  some  pros  and  cons  of  the  various  BCs. 
Section  8  provides  a  short  summary  of  the  BCs. 

Four  appendices  are  included  in  this  report  to  provide  more  detailed  information  on  specific 
aspects  of  BCs  in  the  NCOM.  Appendix  A  provides  a  useful  list  of  the  many  Fortran  variable  names 
used  throughout  the  report.  Appendix  B  presents  details  on  the  NCOM  grid  indexing  scheme  for 
open  boundaries.  Appendix  C  gives  the  advective  BCs  in  finite  difference  form  for  each  of  the 
external  boundaries  (west,  east,  south,  and  north)  for  ease  of  comparison  with  the  Fortran  coding 
within  the  NCOM.  Finally,  Appendix  D  provides  a  glossary  of  acronyms  used  in  this  report. 

2.  NOTATION 

The  lateral  grid  dimensions  in  the  NCOM  are  denoted  by  (n,  m)  where  the  first  index  refers 
to  the  x-direction  and  the  second  index  to  the  y-direction.  Domain  decomposition  has  been  im¬ 
plemented  in  the  NCOM  to  make  it  computationally  scalable  on  parallel  computer  architectures. 
This  means  the  NCOM  grid  can  be  decomposed  into  subsections  (tiles)  so  that  computations  on 
each  subsection  can  be  performed  on  a  separate  processor.  Since  each  tile  requires  information 
from  its  adjacent  neighbors  to  perform  the  finite  difference  computations,  the  implementation  of 
domain  decomposition  requires  that  the  tiles  overlap  each  other  by  one  or  more  grid  points,  i.e. , 
have  a  boundary  zone  or  halo  one  or  more  grid  points  wide  around  each  tile.  The  method  of  domain 
decomposition  is  further  described  in  Wallcraft  and  Moore  (1997). 

Throughout  the  remainder  of  this  report  we  use  AU  to  represent  the  time  step  for  the  model 
(note  that  the  NCOM  currently  uses  the  same  timestep  for  the  internal  and  external  modes).  The 
horizontal  grid  spacing  at  the  mass  points  of  the  C-grid  is  denoted  by  Ax  and  Ay,  and  super¬ 
scripts  (e.g.,  Ax^,  Ax'^)  arc  used  to  denote  the  corresponding  grid  spacing  at  the  velocity  points 
(c.f.  Appendix  B).  D  is  the  total  ocean  depth  (D  =  C  -  H),  where  C  is  the  free  surface  deviation 
and  z  =  H  is  the  depth  of  the  ocean  bottom.  Superscripts  (£)^,D^)  denote  the  corresponding 
values  at  the  velocity  points.  Finite  differences  in  variables  X  are  denoted  by  62tX^  SxX,  and  SyX 
for  the  2At  time  steps,  A.x  and  Ay  spacing,  respectively. 

To  help  the  reader  determine  how  to  specify  the  various  BCs,  we  frequently  specify  the  corre¬ 
sponding  Fortran  variable  name  used  in  the  NCOM.  This  is  given  in  capital  letters,  usually  within 
parentheses  or  brackets  following  the  description  of  the  associated  quantity. 

3.  CYCLIC  BC 

Cyclic  boundaries  in  the  horizontal  x  and  y  directions  are  very  useful  for  running  test  problems 
and  performing  process  studies  to  investigate  the  response  of  simple  ocean  domains  to  idealized 
forcing.  A  common  test  problem  is  a  periodic  alongshore  channel,  which  can  be  used  to  investigate 
processes  occur ing  along  an  idealized  section  of  coast. 
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The  cyclic  boundaries  in  the  NCOM  can  also  be  used  to  perform  1-D  simulations  in  the  vertical 
or  2-D  simulations  within  a  vertical  plane  since  the  cyclic  boundaries  can  maintain  horizontal 
homogeneity  in  the  direction  of  periodicity.  Of  course,  in  this  mode  of  operation  the  NCOM  will 
not  be  as  efficient  as  a  true  1-D  or  2-D  model;  however,  efficiency  is  usually  not  an  important 
factor  in  1-D  and  2-D  studies.  One-dimensional  simulations  can  be  used  to  investigate  surface  or 
bottom  mixed  layer  dynamics  and  the  response  of  coupled  bio-physical  systems  to  vertical  mixing 
and  various  surface  forcings.  Two-dimensional  simulations  can  be  used  to  look  at,  for  example,  the 
dynamics  along  a  single  section  normal  to  a  coast. 

Cyclic  boundaries  in  the  lateral  directions  are  specified  via  the  INDCYC  parameter 
(c.f.  Appendix  A).  The  allowed  possible  values  are 

INDCYC  =  0  (default)  which  specifies  no  cyclic  boundaries, 

INDCYC  =  4  for  cyclic  boundaries  in  the  a:-direction, 

INDCYC  =  5  for  cyclic  boundaries  in  the  y-direction, 

INDCYC  =  6  for  cyclic  boundaries  in  both  the  x  and  y  directions. 

When  performing  1-D  or  2-D  simulations,  the  grid  dimension  in  the  cyclic  direction  is  set  to  one  if 
the  model  is  being  run  with  low  (i.e. ,  2nd-order)  differencing,  but  must  be  set  to  two  if  any  higher- 
order  differencing  is  being  used.  The  grid  overlap  required  for  performing  model  calculations  in  the 
cyclic  direction  is  handled  by  the  halos. 

Note  that  previously  a  provision  existed  within  the  NCOM  for  applying  cyclic  boundaries 
without  using  halos.  In  this  implementation  of  cyclic  boundaries,  three  points  at  the  beginning  of 
the  grid  in  the  cyclic  direction  (points  1  to  3)  were  overlapped  with  three  points  at  the  end  of  the 
grid  (e.g.,  points  n-2  to  n  in  the  x-direction),  and  the  minimum  required  number  of  grid  points 
in  the  cyclic  direction  was  four  (to  allow  a  3-point  overlap).  The  use  of  these  cyclic  boundaries 
was  specified  with  INDCYC  values  from  1  to  3.  However,  the  recent  implementation  of  options  for 
higher-order  differences  within  the  NCOM  requires  a  larger  overlap  in  the  cyclic  grid  direction  (5 
points  rather  than  3).  Hence,  cyclic  boundaries  are  currently  only  implemented  via  halos. 

When  cyclic  boundaries  are  used,  no  open  BCs  need  to  be  specified  in  the  cyclic  direction. 
Note  that  the  halos  are  always  updated  as  if  cyclic  boundaries  are  being  used,  even  if  they  are  not. 
In  the  case  where  closed  or  open  lateral  boundaries  are  being  used,  the  boundary  calculations  are 
not  affected  by  the  halo  values  since  the  latter  are  not  used  in  the  calculations. 

4.  AIR-SEA  BC 

A  variety  of  air-sea  BC  options  are  available  in  the  NCOM  to  provide  the  user  with  some 
flexibility  in  specifying  surface  atmospheric  forcing.  Surface  atmospheric  forcing  fields  that  can  be 
specified  include  surface  air  pressure,  surface  wind  stress,  penetrating  solar  radiation,  and  surface 
fluxes  of  scalar  variables  such  as  heat  and  salt  fluxes. 

OGCMs  frequently  have  difficulty  providing  accurate  simulations  of  sea  surface  temperature 
(SST)  and  sea  surface  salinity  (SSS)  when  driven  only  by  externally  prescribed  fluxes  since  sys¬ 
tematic  errors  in  the  fluxes  or  the  OGCM  will  tend  to  cause  the  model  SST  and  SSS  to  drift  away 
from  observed  values.  For  this  reason,  the  NCOM  can  also  be  forced  by  heat  and  freshwater  fluxes 
calculated  from  a  relaxation  of  the  model  SST  and  SSS  to  externally  prescribed  SST  and  SSS  fields. 
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It  is  generally  preferable  to  relax  the  model  SST  and  SSS  values  to  prescribed  values  by  calcu¬ 
lating  a  flux  rather  than  by  direct  relaxation  to  the  prescribed  values  since  the  relaxation  rate  in 
the  latter  case  depends  on  the  upper-layer  thickness,  which  can  vary  significantly  (e.g.,  over  several 
orders  of  magnitude  if  sigma  coordinates  are  being  used).  If  the  relaxation  is  done  via  a  corrective 
flux,  the  rate  of  relaxation  will  be  independent  of  the  upper-layer  thickness.  Note  that  the  two 
methods  are  formally  equivalent. 

The  flux  that  results  from  relaxation  to  prescribed  values  can  be  used  as  the  sole  flux  or  can  be 
added  to  the  externally  prescribed  flux.  The  latter  procedure  is  preferable  since  with  this  procedure 
the  SST  and  SSS  relaxation  carries  less  of  the  burden  of  determining  the  entire  surface  flux.  For 
the  same  reason,  it  is  important  that  the  externally  prescribed  fluxes  be  accurate,  e.g.,  if  the  model 
and  the  prescribed  fluxes  were  perfect,  no  correction  of  the  model  SST  and  SSS  would  be  needed. 
It  is  especially  important  to  provide  fairly  accurate  estimates  of  the  solar  flux  since  a  significant 
part  of  the  solar  flux  penetrates  below  the  ocean  surface  (which  is  why  it  is  treated  separately  from 
the  surface  heat  flux).  The  penetration  of  solar  radiation  can  have  a  significant  effect  on  the  SST 
(Rochford  et  al.  2001),  upper-ocean  mixing,  and  the  formation  of  the  seasonal  thermocline  (Martin 
1985). 

The  surface  forcing  desired  is  activated  by  setting  the  appropriate  surface  forcing  options  and 
parameters  (c.f.  Appendix  A)  and  providing  the  appropriate  surface  fields.  When  the  input  flag 
INDSBC  =  0,  there  is  no  surface  atmospheric  forcing  and  no  surface  atmospheric  forcing  fields  need 
be  supplied.  This  flag  overrides  the  other  surface  atmospheric  forcing  flags.  When  INDSBC  ==  1, 
surface  atmospheric  forcing  is  applied  and  the  various  surface  atmospheric  fields  can  be  individually 
selected  via  input  flags. 

Most  of  the  flags  for  the  individual  surface  forcing  fields  have  three  possible  values:  0  if  the 
forcing  field  is  not  applied  to  the  ocean,  1  if  the  field  is  read  from  an  input  file,  and  2  if  the  field  is 
to  be  obtained  from  a  coupled  atmospheric  model  (this  last  option  is  not  yet  implemented).  The 
input  flags  for  the  individual  atmospheric  forcing  fields  are  atmospheric  pressure  P  (INDATP), 
wind  stress  r  =  {^xiTy)  (INDTAU),  surface  fluxes  for  temperature  T  (INDSFT)  and  salinity  S 
(INDSFS),  and  surface  solar  irradiance  IR  (INDSOL).  Note  that  each  of  these  surface  forcing 
fields  must  be  provided  to  the  NCOM  on  an  A-grid  configuration. 

There  is  an  additional  input  flag  for  the  solar  radiation  field  (INDCLD)  to  denote  whether  the 
input  field  is  actually  solar  radiation  or  some  representation  of  the  amount  of  cloud  cover.  The 
option  for  cloud  cover  input  is  not  currently  implemented. 

The  NCOM  currently  assumes  that  the  surface  albedo  is  accounted  for  in  the  input  solar 
radiation,  i.e. ,  the  input  solar  radiation  is  the  net  solar  radiation  absorbed  by  the  ocean  and  does 
not  include  the  part  of  the  solar  radation  that  is  reflected  back  into  the  atmosphere. 

The  depth  of  penetration  of  solar  radiation  into  the  ocean  is  governed  by  two  input  flags: 
INDEXT  provides  for  penetration  of  solar  radiation  if  set  to  1  and  does  not  if  set  to  0,  and 
INDTYPE  is  used  to  select  a  solar  extinction  profile.  Values  of  INDTYPE  from  1  to  5  correspond 
to  Jerlov  seawater  optical  types  I,  lA,  IB,  II,  and  III  (Jerlov  1968,  Mobley  1994).  These  optical 
types  range  from  fairly  clear  water  found  in  the  the  central  ocean  gyres  (type  I)  to  fairly  turbid 
water  (type  III).  Note  that  these  input  parameters  allow  only  a  single  seawater  turbidity  type 
within  the  domain.  However,  solar  extinction  within  the  NCOM  is  stored  at  the  upper  interface  of 
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each  grid  cell  in  a  3-D  array  (EXT),  which  is  defined  in  subroutine  SOLEXT.  Hence,  spatially  and 
time  varying  solar  extinction  could  easily  be  implemented  (Rochford  et  al.  2001). 

The  surface  wind  stress,  solar  radiation,  and  scalar  (heat  and  salt)  surface  fluxes  are  defined 
as  positive  downward.  That  is,  a  positive  value  of  r  indicates  a  surface  stress  acting  to  drive  the 
surface  current  in  the  positive  a:  or  y  direction,  and  a  positive  value  of  the  solar  or  surface  heat  flux 
will  act  to  warm  the  surface  layer  of  the  ocean.  Note  that  this  is  the  reverse  of  the  convention  used 
by  the  Princeton  Ocean  Model  (POM),  where  the  surface  heat  fluxes  are  defined  to  be  positive 
upward. 

All  surface  forcing  fields  must  be  supplied  in  dynamical  units,  i.e. ,  the  units  of  the  associated 
physical  quantity  being  forced.  The  surface  atmospheric  pressure  is  expressed  in  terms  of  meters 
of  water.  Note  that  only  the  horizontal  gradient  of  P  drives  the  ocean,  i.e. ,  the  mean  value  of  P 
has  no  effect.  For  a  surface  pressure  given  in  mb,  it  is  suggested  that  the  pressure  in  dynamical 
units  P'  be  calculated  as 


P'  =  {P-  1000)  X  100/(5po),  (1) 

where  g  is  the  gravitational  acceleration  and  po  the  density  of  seawater.  Note  that  an  atmospheric 
pressure  of  1  mb  =  100  Nm'^  =  lOOkgm”^  s“^  corresponds  to  a  surface  elevation  of  about  1  cm. 

For  wind  stress,  the  units  are  m^s“^.  This  is  obtained  from  the  standard  definition  of  wind 
stress  via 

=  r/po,  (2) 

where  the  quantity  u*  is  referred  to  as  the  surface  friction  velocity. 

The  dynamical  units  of  the  surface  fluxes  for  scalar  fields  (i?*tt*)  are  (scalar  units)  ms“^.  For 
the  particular  case  of  surface  heat  and  salt  fluxes,  they  are  °Cms~^  and  psurns”^,  respectively. 
The  temperature  flux  (T*u*)  is  obtained  from 

=  Q/(poC'p),  (3) 

where  Q  is  the  net  surface  heat  flux  (the  sum  of  the  shortwave  and  longwave  radiation  and  latent 

and  sensible  heat  fluxes)  and  Cp  is  the  specific  heat  capacity  of  seawater.  The  salinity  flux  {S,fU^) 
is  obtained  from  the  excess  moisture  flux  of  evaporation  over  precipitation  {E  —  P)  multiplied  by 
the  surface  salinity  5o 


=  5o(E  -  P)/po.  (4) 

The  temporal  variability  of  the  externally  provided  surface  forcing  fields  must  be  the  same  for 
P',  and  This  requires  time  interpolation  or  frequency  subsampling  when  using  fields 

of  different  temporal  variability  within  this  group.  Prescribed  SST  and  SSS  values  can  be  input 
with  different  temporal  variability.  For  real-time  forcing,  the  elapsed  time  of  the  model  run  is  used 
internally  within  the  NCOM  for  the  time  interpolation  of  the  surface  field  to  the  model  time  step. 
For  climate  forcing,  the  time  is  computed  relative  to  the  climate  period  (CLIMATP  =  365  days). 
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The  option  to  relax  the  model  SST  and  SSS  to  specified  values  {Ts  and  Ss^  respectively)  is 
selected  by  setting  the  corresponding  input  flag  (INDSST  or  INDSSS)  equal  to  1  if  the  SST  and 
SSS  are  input  with  the  same  temporal  variability  If  SSS  is  input  with  a  temporal  variability 
different  from  that  used  for  SST,  INDSSS  is  set  equal  to  2.  The  corrective  surface  fluxes  for  the 
relaxation  are  calculated  as 


=  Ar(2o  —  Ts),  (5) 

=  \s{So  -  Ss),  (6) 

where  \t  (RLAXSST)  and  \s  (RLAXSSS)  are  the  rate  of  relaxation  of  the  model  SST  and  SSS  in 
meters  per  day  (md“^),  respectively.  If  prescribed  fluxes  and  relaxation  are  both  being  used,  the 
two  sets  of  fluxes  are  added  together. 

The  surface  atmospheric  forcing  fields  are  defined  within  subroutine  OSURFBC  and  are  linearly 
interpolated  in  time  from  the  prescribed  values. 

5.  SURFACE  ROUGHNESS  AND  TKE 

The  surface  roughness  and  the  surface  value  or  flux  of  turbulent  kinetic  energy  (TKE)  can  be 
specified  for  the  turbulence  submodels.  The  surface  roughness  in  the  NCOM  is  the  effective  vertical 
mixing  length  at  the  surface  of  the  ocean  due  to  wave  action.  The  surface  roughness  is  currently 
used  in  the  turbulence  models  to  provide  a  surface  BC  for  the  vertical  turbulent  mixing  length. 

The  surface  roughness  is  stored  in  a  2-D  array  and  can  be  specified  at  each  horizontal  model 
grid  point.  Integer  flag  INDSURF  is  used  to  specify  the  method  to  be  used  to  determine  the 
surface  roughess.  For  INDSURF  =  0  the  surface  roughness  is  set  to  zero  and  for  INDSURF  = 
2  the  surface  roughness  Z{)  is  estimated  from  the  surface  wind  stress  (friction  velocity)  using  the 
Charnock  relation 


^0  =  QiuVff,  (7) 

where  c/^,  is  the  Charnock  constant  (CHARNOK).  The  value  INDSURF  =  1  is  to  provide  for  2-D 
fields  of  surface  roughness  to  be  read  from  an  input  file,  the  idea  being  to  use  data  from  a  wave 
model  or  a  wave  height  analysis.  This  option  is  not  currently  provided,  but  could  be  implemented 
by  a  user. 

The  Mellor-Yamada  Level  2.5  turbulence  model  (Mellor  and  Yamada  1974,  1982)  also  requires 
a  surface  value  of  TKE  or  a  value  of  the  surface  TKE  flux  as  discussed  by  Craig  and  Banner  (1994). 
Integer  flag  INDTKES  is  used  to  select  which  BC  is  used.  For  INDTKES  =  1,  the  surface  value  of 
TKE  is  calculated  from  the  surface  friction  velocity  as  where  bi  is  a  constant  used  to  scale 

the  dissipation  of  TKE.  For  INDTKES  =  2,  the  surface  flux  of  TKE  is  estimated  from  the  friction 
velocity  as  lOOu^. 

6.  OPEN  BC 

Most  applications  of  an  OGCM  to  an  ocean  basin  or  a  coastal  region  require  that  the  OGCM 
use  information  along  the  open  boundaries  of  the  model  domain  to  achieve  a  realistic  simulation  of 
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the  local  ocean  dynamics.  This  information  may  be  provided  by  larger-scale  OGCMs,  climatologies, 
and  larger  grids  within  the  OGCM  if  it  is  nested. 

The  NCOM  in  its  current  configuration  can  be  run  using  open  boundary  condition  data  from 
the  following  sources:  no  open  boundaries  (INDOBC  =  0),  i.e. ,  a  rigid  wall  along  the  exterior 
boundary  of  the  grid;  initial  values  used  at  boundary  points  for  all  time  (INDOBC  ==  1),  i.e. ,  static 
open  BCs;  boundary  data  from  an  input  file  (INDOBC  =  2)  with  data  either  temporally  varying 
or  static;  and  boundary  data  from  the  coarser  grid  in  which  it  is  nested  in  the  case  of  a  nested 
grid  (INDOBC  =  3).  Note  that  cyclic  boundaries  are  not  included  in  the  grid  points  considered  to 
define  the  open  boundaries. 

It  is  important  for  an  OGCM  to  be  able  to  accommodate  external  information  at  its  open 
boundaries  such  as  surface  elevation,  velocities  and  transports,  values  of  scalar  fields,  and  tides.  The 
velocities  and  transports  at  the  boundaries  must  be  consistent  with  the  density  field.  The  OGCM 
must  also  be  able  to  radiate  interior  waves  out  of  the  model  domain  with  minimal  consequences 
on  the  OGCM  simulation.  A  variety  of  open  BCs  have  been  implemented  in  the  NCOM  based  on 
findings  from  other  modelers  that  have  demonstrated  good  performance  in  fulfilling  these  criteria. 

Open  BCs  are  specified  via  the  OPENBC  subroutine  of  the  NCOM,  Boundary  conditions  have 
been  implemented  for  the  elevation,  barotropic  normal  and  tangential  velocity,  baroclinic  normal 
and  tangential  velocity,  vertical  velocity,  scalar  fields  (e.g.,  temperature  and  salinity),  and  vertical 
mixing  variables.  These  are  each  described  in  turn  below.  Note  that  some  of  the  BCs  work  better 
than  others.  Some  pros  and  cons  of  the  BC  are  discussed  in  this  section  and  in  Section  7.  The 
Fortran  variable  names  for  the  various  fields  can  be  found  in  Appendix  A, 

6.1  Elevation  and  Normal  Barotropic  Velocity 

Two  BCs  are  provided  for  dealing  with  the  surface  elevation  and  normal  barotropic  velocity 
at  the  open  boundaries:  clamped  (INDOBE  1)  and  Flather  (INDOBE  =  2).  These  BCs  involve 
both  the  surface  elevation  and  the  normal  barotropic  velocity  at  the  boundary;  hence,  the  BC  for 
these  two  fields  are  discussed  together. 

6.1.1  Calculation  of  the  Free- Surface  Mode 

To  understand  how  the  elevation  and  barotropic  BCs  influence  the  model  solution,  a  brief 
explanation  of  the  free-surface  mode  calculation  is  warranted.  In  the  NCOM,  the  free  surface  is 
calculated  implicitly  where  the  surface  pressure  gradients  in  the  momentum  equations  and  the 
divergence  terms  in  the  surface  elevation  equation  have  a  component  at  the  new  time  level  (n  + 1) 
being  calculated.  The  finite  difference  equations  for  the  free-surface  mode  are  (Martin  2000) 

+  asC”-')  +  (8) 

Ax^D'^gSyiaiC^^  +  aaC"  +  (9) 

<5,(A2/“(/?i(D"u)"+i  +  f32{D^ur  + 

-  6y{Ax^{pi{D^vr+^  +  p2{D^vr  +  03{D^vr-^))  +  AxAyDQ,  (10) 


Ax^Ay'" 

2Ai 

Ax^Aif 


2At 


S2i{D'‘u)  =  - 
52t{D^v)  =  - 


AxAy 

2At 


52iC  =  - 
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where  (u^v)  are  the  depth  averaged  velocities.  The  D^Gu  and  D^Gy  are  the  vertical  integrals  of 
all  the  forcing  terms  of  the  equations  for  the  baroclinic  velocity  u  and  u,  respectively,  except  for 
the  surface  elevation  gradient  terms,  and  and  D'^  —  .  Q  is  the  depth-averaged  mass 

flux  source  term.  The  ai  (Pi)  specify  the  fractional  weighting  of  the  surface  elevation  gradient  in 
the  momentum  equations  (divergence  terms  in  the  depth-averaged  continuity  equation)  at  the  new 
(i  =  1),  current  (i  —  2),  and  previous  (z  =  3)  time  levels.  See  Section  3.6  of  Martin  (2000)  for 
further  details. 

The  equations  for  the  free-surface  mode  are  solved  by  substituting  the  expressions  for 
and  from  Eqs.  (8)  and  (9)  into  Eq.  (10)  and  solving  for  the  new  surface  elevation 

The  resulting  equation  is  an  elliptic  equation  for  which  can  be  solved  with  an  iterative  or  a 

direct  method  (NCOM  currently  uses  an  iterative  solver). 

The  free  surface  mode  is  solved  such  that  the  baroclinic  velocities,  barotropic  transports,  and 
elevation  are  self-consistent  with  respect  to  each  other.  This  is  accomplished  via  the  following 
calculation  sequence: 

1.  Advection  velocities  and  horizontal  eddy  coefficients  needed  for  momentum,  new  densities,  and 
new  baroclinic  pressure  gradients  are  calculated. 

2.  New  3-D  horizontal  velocities  are  calculated  and  the  forcing  terms  from  the  3-D  momentum 
equations  are  vertically  integrated  to  provide  the  forcing  terms  needed  for  the  depth-averaged 
momentum  equations  that  are  used  to  calculate  the  free-surface  mode. 

3.  The  depth-averaged  momentum  and  continuity  equations  are  solved  for  the  new  surface  elevation 
and  depth- averaged  velocities. 

4.  The  new  3-D  velocity  fields  calculated  in  Step  2  are  corrected  by  adding  a  depth-independent 
adjustment,  so  that  their  vertical  mean  agrees  with  the  new  depth- averaged  velocities  calculated 
in  Step  3.  This  effectively  corrects  the  3-D  velocities  for  the  new  surface  elevation  gradient. 

5.  The  velocity  field  that  will  be  used  to  advect  the  scalar  fields  is  calculated  by  adding  a  depth- 
independent  adjustment  to  the  3-D  velocity  fields  at  time  level  n,  so  that  the  depth- averaged 
advection  velocities  are  consistent  with  the  depth- averaged  continuity  equation. 

6.  An  Asselin  filter  is  applied  to  the  velocity  and  surface  elevation  fields.  The  filtered  3-D  velocities 
are  then  corrected  to  be  consistent  with  the  filtered  depth-averaged  velocities  using  the  same 
procedure  as  in  Step  4. 

The  adjustment  of  the  advection  velocities  in  Step  5  ensures  that  the  velocity  field  used  to 
advect  the  scalar  fields  is  numerically  nondivergent.  This  is  necessary  to  avoid  spurious  sources 
and  sinks  when  using  the  flux  form  of  numerical  advection. 

6.1,2  Clamped  BC 

With  a  clamped  BC  (INDOBE  =  1)  the  boundary  condition  used  for  the  free  surface  mode  of 
the  model  is  the  specified  surface  elevation  at  the  boundary.  This  BC  is  generally  used  only  when 
tidal  forcing  is  the  sole  forcing  both  at  the  boundary  and  within  the  model  domain.  The  reason 
for  this  restriction  is  that  surface  waves  generated  by  other  forcing  within  the  model  domain  will 
reflect  from  the  boundary  back  into  the  interior  when  a  clamped  BC  is  used.  These  reflections  can 
severely  or  catastrophically  degrade  the  model  solution  in  the  interior  of  the  domain. 
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Transient  surface  waves  generated  as  the  tidal  cycle  spins  up  will  also  reflect  at  the  boundary 
These  transient  waves  will  eventually  dissipate  due  to  internal  damping,  mainly  from  bottom  fric¬ 
tion.  The  time  required  for  the  transient  waves  to  die  away  and  for  the  model  solution  to  settle 
down  to  a  repeating  tidal  cycle  will  depend  on  the  bottom  depth  and  the  magnitude  of  the  bottom 
friction,  and  can  range  from  a  day  or  two  for  a  shallow  coastal  region  to  a  couple  of  months  for  a 
large  basin. 

With  a  clamped  BC  the  normal  baroclinic  velocity  at  the  boundary  is  calculated  by  the  model; 
hence,  data  for  this  field  at  the  boundary  are  not  needed.  Note  that  tidal  simulations  with  a 
clamped  BC  can  provide  a  means  for  calculating  the  tidal  velocities  at  the  boundary  that  are 
needed  for  using  the  Flather  radiation  BC. 

6.1.3  Flather  BC 

With  the  Flather  BC  (INDOBE  =  2),  the  BC  is  on  the  normal  barotropic  transport  (14  =  VnD, 
where  Vn  is  the  depth-averaged  normal  velocity)  rather  than  the  elevation  (Flather  and  Proctor 
1983).  The  advantage  of  the  Flather  BC  over  the  clamped  BC  is  that  the  Flather  BC  allows 
barotropic  disturbances  within  the  model  domain  to  propagate  out  through  the  boundaries  (i.e. , 
radiate)  rather  than  reflect  back  into  the  interior.  The  disadvantage  is  that  values  of  both  the 
surface  elevation  and  the  normal  barotropic  transport  are  needed  for  implementation.  The  normal 
barotropic  transport  for  the  Flather  BC  is  calculated  within  the  NCOM  as 

v;,  =  +  Cbtiv  -  (n) 

where  and  are  the  prescribed  normal  barotropic  transport  and  elevation  at  the  boundary 
for  the  ocean  circulation,  and  are  the  normal  barotropic  transport  and  elevation  for  the 
prescribed  tidal  forcing,  Cu  is  the  surface  wave  speed 

Cu  =  (12) 

and  rj  is  the  NCOM  elevation  at  the  first  interior  point  next  to  the  boundary. 

The  Flather  BC  is  quite  flexible  because  it  allows  both  tidal  and  circulation  transports  to  be 
specified  at  the  boundary,  while  at  the  same  time  allowing  barotropic  disturbances  to  propagate 
out.  It  can  sometimes  be  used  to  relax  and  to  interior  values  in  the  case  where  values  of 
and  77“’’^  are  unavailable.  We  note  that  the  Flather  radiation  condition  has  also  been  applied 
in  other  OGCMs  (Oey  and  Chen  1992a-b;  Rochford  and  Shulman  2000). 

6.2  Tangential  Barotropic  Velocity 

The  BCs  for  the  tangential  barotropic  velocity  are  formulated  in  terms  of  the  tangential 
barotropic  transport  (V  =  vtD  where  vt  is  the  depth-averaged  tangential  velocity).  There  are 
three  BCs  provided  for  the  tangential  barotropic  velocity;  a  zero  gradient  (INDOBVB  =  1),  an 
Orlanski  radiation  condition  (INDOBVB  =  2),  and  an  advective  BC  (INDOBVB  =  3). 

6.2.1  Zero  Gradient  BC 

When  a  zero  gradient  condition  is  applied  (INDOBVB  =  1)  the  tangential  barotropic  transport 
at  the  open  boundary  is  set  equal  to  the  value  at  the  first  interior  point  next  to  the  boundary,  i.e. , 
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the  field  is  set  to  have  a  zero  gradient  in  the  direction  normal  to  the  boundary.  This  is  a  simple 
BC,  but  it  frequently  works  fairly  well. 

6.2.2  Orlanski  Radiation  BC 


For  the  Orlanski  radiation  BC  option  (INDOB VB  =  2),  the  NCOM  uses  an  Orlanski  radiation 
BC  for  outward  propagating  barotropic  disturbances  and  relaxation  to  an  externally  prescribed 
value  for  inflowing  disturbances.  The  radiation  condition  introduced  by  Orlanski  (1976)  (see  also 
Miller  and  Thorpe  (1981))  is  a  variation  of  the  Sommerfeld  radiation  condition.  For  the  tangential 
barotropic  transport  the  equation  for  the  Orlanski  radiation  condition  is 


dt  ^ dXn 


=  0, 


(13) 


where  the  phase  speed  c  is  calculated  for  the  variable  Vt  at  each  boundary  point  and  is  the 
coordinate  normal  to  the  boundary  with  increasing  value  in  the  inflow  direction.  Note  that  the 
phase  speed  is  deflned  to  be  negative  for  inflow  (c  <  0)  and  positive  for  outflow  (c  >  0) . 


The  Orlanski  radiation  condition  is  implemented  in  the  NCOM  using  a  leapfrog  temporal 
formulation  as 


Vt{B,  N  +  \)^  [(1  -  a)  Vt{B,  N-l)  +  2a  Vt{B  -  1,  iV)] ,  (14) 

where  B  denotes  the  grid  point  location  of  the  boundary,  B  ~1  denotes  the  first  interior  grid  point 
location  from  the  boundary,  N  denotes  the  time  level,  and  a  is  given  by 

a  =  max[0,  min[l,  cAt/Ax]],  (15) 


where 


cAt  ^  Vt{B  -  1,  TV  -  1)  -  V,{B  -  1,  TV  +  1) 

Ax  Vi{B-l,N-\-l)-^  Vt{B  2Vt{B  ^  ^ 

and  requires  values  of  Vt  at  the  three  time  levels  N  ~  and  N  ~^1. 

The  Orlanski  radiation  BC  is  only  used  for  an  outward  propagating  signal  (c  >  0).  The 
relaxation  in  the  case  of  an  inward  propagating  signal  (c  <  0)  is  governed  by 

-^  =  Xv  {vr^  +  -  Vt)  ,  (17) 

where  Vt^^^  is  the  prescribed  tangential  barotropic  velocity  for  the  ocean  circulation,  is  the  pre¬ 
scribed  tidal  contribution,  and  is  the  time  scale  of  relaxation  for  Vt  [RLXOBVB  =  (AyAti)"^]. 

6.2.3  Advective  BC 

An  alternative  approach  to  the  Orlanski  radiation  BC  is  to  apply  one-way  external  forcing  using 
an  advective  BC  (INDOBVB  —  3)  for  the  tangential  barotropic  transports  (Rochford  and  Shulman 
2000).  This  is  a  forward,  upwind,  two-time-level  scheme  where  the  advection  term  is  lagged  (i.e. ,  is 
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at  the  central  time  level).  Here  the  externally  specified  tangential  barotropic  transport  is  advected 
into  the  model  domain  in  the  case  of  inflow,  and  the  internal  tangential  barotropic  transport  is 
advected  to  the  open  boundary  in  the  case  of  outflow.  Note  that  we  define  Un  ^  0  for  inflow  and 
<  0  for  outflow  where  Vn  is  the  depth-averaged  normal  velocity.  This  leads  to  a  BC  of  the  form 


dVt 

dt 


dVt 

'dZyi 


,  -  A 


(18) 


where 

dVt  ^  f  {Vt-  Vob)I^Xn  >  0  (inflow) 
dXn  ~  \  {Vi-  Vt)I^Xn  Vn  <  0  (outflow), 


(19) 


with  Vi  the  tangential  barotropic  transport  at  the  first  grid  point  inside  the  open  boundary  \yt{B  — 
1,N)],  Vob  the  tangential  barotropic  transport  provided  at  the  open  boundary,  and  Ax„  the  grid 
spacing  in  the  direction  normal  to  the  boundary. 


The  advective  BC  is  implemented  in  the  NCOM  using  a  leapfrog  temporal  formulation 


Vt{B,  N  +  l)  =  Vt{B,  N)  -  +  \c\)[Vt{B,  N)  -  Vob{B,  iV)]  + 

"  (c-|c|)[Hi(B,iV)-V^(B,iV)]}, 

VniB,N) 
max.{€,  D{B,N))' 


(20) 


(21) 


where  AU  is  the  model  time  step  and  e  =  10"®  is  included  to  avoid  division  by  zero.  Further  details 
on  Eq.  (19)  can  be  found  in  Appendix  C. 


6.3  Normal  Baroclinic  Velocity 

Four  BCs  are  provided  for  the  normal  baroclinic  velocities;  use  the  value  calculated  internally 
by  the  model  (INDOBU  =  1),  use  an  Orlanski  radiation  condition  (INDOBU  =  2),  combine  an 
upwind  horizontal  advection  to  the  value  calculated  by  the  model  (INDOBU  =  3),  and  use  an 
advective-radiative  BC  (INDOBU  =  4). 

Note  that  all  baroclinic  model  velocities,  including  the  values  at  the  boundaries,  are  adusted 
at  each  timestep  at  each  horizontal  point  by  a  constant  value  such  that  the  depth-averaged  value 
agrees  with  the  depth-averaged  velocity  computed  for  the  barotropic  mode  (c.f.  Section  6.1.1).  This 
means  barotropic  signals  (e.g.,  tidal  signals)  will  be  incorporated  into  the  baroclinic  velocities  at  the 
open  boundary  points  even  though  they  are  not  explicitly  accounted  for  when  setting  the  baroclinic 
velocities  at  the  boundary  points. 

6.3.1  Internal  Model  Calculated  Value  for  BC 

The  NCOM  always  calculates  a  value  for  the  normal  baroclinic  velocity  at  the  open  boundary 
points  using  the  same  calculation  as  for  the  interior  velocity  points,  except  that  the  advection 
and  horizontal  mixing  terms  are  omitted  (because  all  the  quantities  needed  to  calculate  these 
terms  in  the  manner  of  the  other  interior  points  are  not  available  at  the  open  boundary  points). 


12 


Rockford  and  Martin 


When  INDOBU  ~  1,  this  internally  calculated  value  of  the  normal  baroclinic  velocity  is  taken  as 
the  boundary  value.  Note  that  this  BC  works  fairly  well,  but  is  usually  improved  if  horizontal 
advection  is  included  (i.e. ,  for  INDOBU  =  3). 

6,3,2  Orlanski  Radiation  BC 


For  INDOBU  =  2,  an  Orlanski  radiation  condition  is  used  to  set  the  normal  baroclinic  velocity 
at  the  open  boundary  points  for  outward  propagating  disturbances,  and  a  relaxation  to  an  externally 
prescribed  value  is  used  for  inward  propagating  disturbances.  These  BCs  are  calculated  in  similar 
fashion  to  those  for  the  tangential  barotropic  velocity.  The  time  scale  for  relaxation  to  the  externally 
prescribed  normal  baroclinic  velocity  is  specified  by  the  input  parameter  RLXOBV, 

6.3.3  Internal  Model  Calculated  Velocity  plus  Advection  BC 


For  INDOBU  =  3,  the  normal  baroclinic  velocity  calculated  internally  by  the  model  at  the  open 
boundary  points  is  modified  by  adding  to  it  a  change  due  to  advection.  An  upwind  advection  term 
normal  to  the  boundary  is  applied  that  uses  the  interior  model  normal  velocity  when  advection  is 
outward  and  the  prescribed  normal  velocity  when  the  advection  is  inward.  The  equation  describing 
the  advection  is 


dvn  dvn 

dt  dxn 


-0, 


(22) 


where  the  normal  gradient  of  the  velocity  is  calculated  as 

dVn  _  I  (vn  -  Vab)/ALXn  Vn  >  0  (inflow) 
dXn  ~  \  (Vi-  Vn)/AXn  Vn  <  0  (outflow), 


(23) 


where  vi  is  the  baroclinic  velocity  at  one  grid  point  inside  the  open  boundary  [vt{B  —  1,  A”)],  and 
Vob  is  the  specified  (external  or  coarse  grid)  value  at  the  open  boundary. 


This  advective  BC  is  implemented  using  a  leapfrog  scheme  that  is  identical  in  form  to  that 
used  for  the  tangential  barotropic  transports 

N  +  l)=  vtiB,  N)  -  :^{{c  +  \c\)[vt{B,  N)  -  v,,{B,  TV)]  +  (24) 

2iL\Xn 

where  in  this  case  c  ~  Vn{B,N),  and  is  the  internal  time  step.  Further  details  on  Eq.  (23)  are 
found  in  Appendix  C. 


6.3.Jf  Advective-Radiative  BC 


A  radiative  variation  of  the  advective  BC  has  been  found  to  work  quite  well  in  the  POM 
when  applied  to  the  normal  baroclinic  velocity  (Rochford  and  Shulman  2000).  To  compare  NCOM 
simulations  with  those  from  the  POM,  this  same  BC  has  been  included  as  an  option  in  the  NCOM 
(INDOBU  =  4).  The  radiational  BC  is  applied  to  the  normal  baroclinic  velocity  Vn  as 


dvji 


^  -n 

^bc  O  - 

UXji 


(25) 
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where  C^c  is  the  phase  speed  for  outward  radiation  of  the  normal  velocity,  and  the  normal  coordinate 
is  positive  in  the  direction  of  inflow.  While  it  is  possible  to  implement  an  implicit  algorithm 
to  evaluate  the  phase  speed,  a  value  for  C^c  is  explicitly  prescribed  for  computational  efficiency. 
This  requires  that  the  user  specify  C^c  in  units  of  ms"^  via  the  cgwb(ib,2)  variable  in  subroutine 
OPENBC.  We  note  that  the  PWC  POM  of  CoBALT  uses  Cbc  =  0.005  (Rochford  and 

Shulman  2000).  The  optimal  choice  for  Cbc  must  be  determined  through  numerical  tests.  To  reduce 
the  development  of  noise,  a  Hanning  filter  is  applied  afterwards  to  Vn  parallel  to  the  boundary 

Vn{B,N  +  l)  =  {l-Xi)  Vn{B,N)  +  Cbc^^^^^^^  +XlVob  (26) 

Vn{B  -  1,  AT  +  1)  =  (1  -  A2)  Vn{B  -  1,  AT)  +  A2V06,  (27) 

with  Ai  =  0.1  and  A2  =  0.05.  The  normal  derivative  is  given  by  Eq.  (23)  and  acquires  the  same 
finite  difference  form  as  given  in  Eq.  (24). 

6.4  Tangential  Baroclinic  Velocity 

The  open  BCs  available  for  the  tangential  baroclinic  velocities  are  zero  gradient  (INDOBV  =  1), 
Orlanski  radiation  (INDOBV  =  2),  and  advective  (INDOBV  =  3). 

6.4^1  Zero  Gradient  BC 

The  zero  gradient  BC  for  the  tangential  baroclinic  velocity  (INDOBV  =  1)  is  applied  in  similar 
fashion  to  the  zero  gradient  BC  for  the  tangential  barotropic  velocity.  The  values  at  the  boundary 
are  set  equal  to  the  values  of  the  tangential  baroclinic  velocity  just  inside  the  boundary.  This  simple 
BC  frequently  works  fairly  well. 

6,4 Orlanski  Radiation  BC 

The  Orlanski  radiation  BC  for  the  tangential  baroclinic  velocity  (INDOBV  =  2)  is  applied  in 
the  same  way  as  the  Orlanski  radiation  condition  for  the  tangential  barotropic  velocity.  A  discussion 
of  this  BC  will  therefore  not  be  repeated  here. 

5.^.5  Advective  BC 


where  Vi  is  the  baroclinic  velocity  at  one  grid  point  inside  the  open  boundary  —  1,  V)],  and 

Vob  is  the  specified  (external  or  coarse  grid)  value  at  the  open  boundary.  This  BC  is  implemented 
in  the  same  manner  as  for  the  tangential  barotropic  transport  and  normal  baroclinic  velocity 
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-{{c+\c\){vt{B,N)-Voi{B,N)]  + 
'  ic-\c\)[vi{B,N)-vt{B,N)]}, 


with  the  difference  that  here  c  —  Vn{B ^  N)  and  Ati  is  the  internal  time  step. 

6.5  Vertical  Velocity 

A  zero  gradient  BC  in  the  horizontal  is  applied  to  the  vertical  velocity  w  along  all  open  bound¬ 
aries,  i.e. , 


where  Xn  is  in  the  direction  normal  to  the  open  boundary.  With  this  BC,  the  vertical  velocity  at  a 
boundary  point  is  set  equal  to  the  value  at  the  adjacent  interior  grid  point. 

6.6  Scalar  Fields 

Four  BCs  are  provided  for  the  scalar  (e.g.,  temperature  and  salinity)  fields:  the  externally 
prescribed  values  are  put  directly  on  the  boundary  (INDOBR  =1),  Orlanski  radiation  (INDOBR  = 
2),  horizontal  advection  (INDOBR  =  3),  and  full  advection  (INDOBR  ==  4). 

6.6.1  Prescribed  Values 

With  this  BC  (INDOBR  =  1),  the  externally  prescribed  values  of  the  scalar  fields  are  placed 
directly  on  the  boundary.  This  direct  specification  of  the  boundary  values  tends,  in  general,  to  be 
reflective,  i.e. ,  not  radiative.  This  BC  is  used  when  performing  two-way  nesting  when  values  of  the 
scalar  fields  are  fed  back  from  the  nested  grid  to  the  main  grid.  In  this  case,  the  solutions  on  the 
nested  grid  and  on  the  main  grid  are  closely  coupled  and  the  direct  placement  of  the  values  of  the 
scalar  fields  from  the  main  grid  onto  the  boundary  of  the  nest  helps  to  maintain  the  close  coupling. 
Hence,  a  radiative  BC  is  not  required. 

6.6.2  Orlanski  Radiation  BC 

With  this  open  BC  (INDOBR  =  2),  scalars  r  are  relaxed  to  an  externally  prescribed  value  for 
inflow  and  an  Orlanski  radiation  condition  is  applied  for  outflow.  The  relevant  equation  in  the  case 
of  inflow  is 

dr 

~  ^7'  (T oh  ) 


where  rob  is  the  prescribed  scalar  field  at  the  open  boundary,  and  A,.  ^  is  the  corresponding  time 
scale  of  relaxation  [RLXOBR  =  (ArAti)"^].  For  outflow,  the  radiation  condition  applied  is 

dr  dr 

(jt  (JXji 


where  the  phase  speed  c  is  calculated  for  each  scalar  field  r,  and  Xn  is  the  coordinate  normal  to  the 
boundary  with  increasing  value  in  the  inflow  direction. 
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The  relaxation  condition  in  the  case  of  inflow  is  implemented  using  a  leapfrog  temporal  formu¬ 
lation  as 

r{B,  N+l)  =  r{B,  N)  +  2A,  [rob  -  r{B,  A^)] .  (32) 

For  outflow,  the  Orlanski  radiation  condition  for  scalars  is  implemented  as  given  by  Eqs.  (14)  to 
(16)  with  Vt  replaced  by  r.  An  Asselin  temporal  filter  is  applied  after  the  boundary  value  is  updated 
to  prevent  time  splitting  with  the  leapfrog  temporal  integration  scheme. 


6.6.3  Advective  BC 


A  horizontal  advective  BC  is  applied  to  the  scalars  at  the  open  boundaries  that  is  forward 
in  time  (2  time  level)  and  upwind  as  implemented  for  the  barotropic  and  baroclinic  transports 
(INDOBR  =  3) 


dr  dr 

ut  OXfi 


=  0. 


(33) 


The  normal  gradients  of  the  scalars  are  given  here  by 

dr  _  f  (r  -  rob)/^Xn  >  0  (inflow) 
d^  ~  [  In-  r)/Axn  Vn  <  0  (outflow), 


(34) 


where  n  is  the  scalar  at  one  grid  point  inside  the  open  boundary  [r{B  -  1,N)],  and  rob  is  the 
prescribed  (external  or  coarse  grid)  value  at  the  open  boundary.  For  scalars,  this  BC  is  implemented 
in  the  NCOM  as 


r{B,  N)  -  {cAti/Axn)  [r{B,  N)  -  robiB,  N)]  c  >  0  (inflow) 
r(B,  N)  -  (cAti/Axn)  h(B,  N)  -  r{B,  N)]  c  <  0  (outflow). 


(35) 


with  c  =  Vn{B,N).  To  prevent  time  splitting  with  the  leapfrog  temporal  integration  scheme,  an 
Asselin  temporal  filter  is  applied  after  the  boundary  value  is  updated. 


6.6.4  Advective  BC  with  Vertical  Advection 


In  some  situations,  such  as  with  internal  wave  propagation,  or  strong  upwelling  or  downwelling 
along  the  open  boundaries,  it  may  be  desirable  to  include  the  effect  of  vertical  advection  as  well  as 
horizontal  advection  on  the  boundary  values  of  the  scalar  fields.  The  equation  for  vertical  advection 
of  a  scalar  field  at  the  boundary  is 


dr  dr 


0, 


(36) 


where  w  is  the  vertical  velocity  (positive  value  upwards).  Such  an  option  is  available  in  the  NCOM 
(INDOBR  =  4)  as  a  forward  in  time  (2  time  level),  centered  scheme  that  is  applied  after  the  scalars 
have  been  updated  for  the  horizontal  BC  (r') 

r(B,N  +  l,k)  =  r'(B,  N,k)+c  [r(B,  N,k  +  1)-  r(B,  N,k-  1)] 

Ati  w{B,  N,  k)  +  w{B,  N,k+  1) 

2  .  Azfe  +  Azk+i 


(37) 

(38) 
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Here  k  denotes  the  vertical  level  and  =  zw}^  —  zwi^^i  is  the  distance  between  the  top  of  levels 
k  and  k  +  1.  As  in  the  case  with  no  vertical  advection,  an  Asselin  temporal  filter  is  applied  after 
the  boundary  value  is  updated  (c.f.  Section  6.6.3). 

6.7  Vertical  Mixing  Variables 

A  zero  gradient  BC  is  applied  to  the  turbulent  length  scale  (/)  and  the  vertical  eddy  coefficients 
for  momentum  {km)  and  scalars  (kh)  as  in  Section  6.5.  When  the  Mellor-Yamada  level  2.5  scheme 
is  used  (INDZK  =  3),  this  BC  is  also  applied  to  twice  the  turbulent  kinetic  energy  (q^)  and  its 
product  with  the  turbulence  length  scale  L 

6.8  Tidal  Forcing 

The  NCOM  provides  the  option  of  superimposing  tidal  forcing  along  the  open  boundaries  on 
any  boundary  data  that  are  used  (INDTIDE  >1).  This  option  is  typically  used  for  cases  where 
tidal  forcing  is  important  but  is  not  included  in  the  boundary  data,  or  where  some  particular  tidal 
constituents  are  desired  but  are  not  included  in  the  boundary  data.  The  default  option  is  no  tidal 
forcing  beyond  that  implicitly  included  in  the  boundary  fields  (INDTIDE  =  0). 

Tidal  forcing  is  prescribed  at  the  boundary  via  an  input  file  that  specifies  the  particular  tidal 
constituents  (n  =  1, . . .  ^ritc)  for  which  data  are  being  provided  and  the  amplitude  and  phase  of 
the  fields  for  each  of  the  tidal  constituents.  The  fields  presently  included  are  the  elevation  and  the 
normal  and  tangential  barotropic  transport  (velocity  times  depth). 

The  total  elevation  and  barotropic  transports  due  to  the  tides  are  calculated  as  the  sum  of  the 
values  of  the  individual  constituents,  i.e. , 

ntc 

n=l 
ntc 

n~l 
ntc 

VUd.^Y^Vtcos{oji^t-4^), 

n=l 

where  is  the  frequency  of  a  tidal  constituent  in  units  of  Hz  and  t  is  the  elapsed  model  time 
in  seconds.  The  tidal  velocity  transport  variables  are  only  used  if  the  Flather 

radiation  BC  is  applied  for  the  elevation.  They  are  not  used  when  the  elevation  is  clamped  along 
the  open  boundary.  These  latter  two  BCs  are  described  next. 

The  input  tidal  data  should  consist  of  the  equilibrium  tidal  amplitudes  and  phases.  The 
amplitude  and  phase  correction  for  the  individual  tidal  constituents,  which  depends  on  the  location 
and  time  of  year,  will  be  calculated  by  the  NCOM  (in  subroutine  TIDE_FAC)  based  on  the  longitude 
and  latitude  of  the  domain  and  the  initial  date  (IDATE)  and  time  (ITIME)  of  the  simulation. 

Note  that  rotation  of  the  tidal  velocity  data  to  adjust  them  to  the  local  orientation  of  the  model 
grid  is  not  just  a  simple  rotation  since  both  the  amplitudes  and  phases  of  the  velocities  have  io  bo 
taken  into  account. 


(39) 

(40) 

(41) 
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ThG  tides  in  the  ocean  are  forced  by  the  local  tidal  potential,  which  is  due  to  astronomical 
forces.  For  small  coastal  domains  of  less  than  100  to  200  km  in  horizontal  extent,  the  effect  of  the 
local  tidal  potential  on  a  simulation  of  the  tides  is  generally  very  small  and  is  usually  neglected, 
i.e. ,  the  forcing  of  the  tide  at  the  lateral  boundary  of  the  domain  is  sufficient.  However,  the  tidal 
potential  forcing  becomes  increasingly  important  as  the  size  and  depth  of  the  domain  is  increased, 
and  may  need  to  be  accounted  for  in  larger  domains. 

The  tidal  potential  forcing  enters  the  model  in  the  momentum  equations  as  the  horizontal 
derivative  of  the  equilibrium  surface  elevation  corresponding  to  the  tidal  potential  for  a  particular 
tidal  constituent  (this  is  the  same  way  that  the  surface  atmospheric  pressure  is  implemented). 
NCOM  provides  a  variable  (EP)  for  the  total  tidal  potential,  a  subroutine  (TIDEPOT)  to  calculate 
EP,  and  a  logical  flag  (TIDPOT)  to  turn  the  tidal  potential  forcing  on  or  off.  However,  subroutine 
TIDEPOT  is  currently  only  acting  as  a  placeholder  for  the  tidal  potential  calculation,  i.e. ,  it 
currently  just  returns  a  value  of  zero  for  EP.  In  order  to  implement  the  tidal  potential  within  the 
NCOM,  subroutine  TIDEPOT  must  be  modified  to  calculate  the  tidal  potential  corresponding  to 
the  tidal  constituents  being  used.  Note  that  the  tidal  potential  for  a  given  tidal  constituent  is  a 
function  of  the  location  and  time. 

7,  DISCUSSION 

As  noted  previously,  some  BCs  work  better  than  others.  Some  of  the  BCs  in  the  NCOM  have 
been  implemented  for  certain  specific  applications,  or  for  testing  and  comparison  purposes,  and 
may  not  be  good  BCs  for  general  use.  Here  we  provide  some  discussion  of  the  pros  and  cons  of  the 
various  BCs  and  some  recommendations. 

7.1  Air-Sea  BC 

The  NCOM  provides  for  forcing  by  the  surface  atmospheric  pressure  gradient.  This  forcing 
tends  to  be  less  important  than  the  wind  stress  and  fluxes  of  temperature  and  salinity  for  driving 
the  ocean.  For  this  reason  it  is  frequently  not  accounted  for  in  ocean  modeling.  However,  the  effect 
of  atmospheric  pressure  is  not  negligible  and  can  be  important  in  certain  situations,  e.g.,  surface 
pressure  differences  across  straits  can  force  significant  flows  through  the  straits. 

The  surface  wind  stress  is  of  prime  importance  in  driving  the  ocean  and  accurate  wind  stresses 
are  needed  to  get  accurate  predictions  of  ocean  currents  and  mixing.  Obtaining  accurate  wind 
stresses  can  be  a  problem  and  this  is  especially  true  in  coastal  areas  where  available  wind  stress 
data  may  not  be  adequate  to  describe  the  true  temporal  and  spatial  variability. 

For  scalar  fields,  both  surface  fluxes  and  surface  values  can  be  prescribed.  Using  prescribed 
fluxes  alone  tends  to  cause  drift  of  the  model  near-surface  values  over  time  due  to  errors  in  the 
fiuxes  and  in  the  model  itself.  The  drift  can  occur  quite  rapidly,  especially  for  the  temperature  in 
shallow  coastal  areas,  which  can  heat  or  cool  quickly.  Without  feedback  of  the  local  SST  to  the 
surface  heat  flux,  a  prescribed  heat  flux  cannot  respond  to  the  changing  SST  the  way  that  it  should 
and  unrealistic  SST  values  can  develop  within  days  or  even  hours. 

Using  relaxation  of  scalar  fields  to  prescribed  surface  values  alone  may  also  be  undesirable 
since  this  forces  the  relcixation  to  carry  the  full  burden  of  computing  the  surface  flux.  Under  such 
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circumstances  a  short  relaxation  time  scale  will  be  needed  to  keep  the  model  values  close  to  the 
prescribed  values.  A  better  approach  is  to  provide  the  best  estimates  of  the  surface  fluxes  that 
are  available  while  simultaneously  relaxing  to  prescribed  values.  With  this  combination  of  forcing, 
if  the  fluxes  are  reasonably  accurate,  only  a  moderate  rate  of  relaxation  to  the  prescribed  values 
should  be  needed  to  keep  the  model  surface  values  on  track.  It  is  especially  important  to  prescribe 
accurate  solar  radiation  since  solar  heat  penetrates  into  the  sea  and  this  penetration  significantly 
affects  upper  ocean  heating  and  mixing  and  the  development  of  the  seasonal  thermocline. 

An  alternative  for  calculating  surface  fluxes  is  to  use  the  model  SST  and  SSS  along  with  pre¬ 
scribed  surface  marine  parameters  such  as  wind  speed,  cloud  cover,  humidity,  and  air  temperature 
and  pressure  to  calculate  the  fluxes  via  bulk  aerodynamic  formulas  (Rochford  et  al.  2000,  Kara  et 
al.  2000).  The  surface  marine  parameters  can  be  obtained  from  atmospheric  model  forecasts  in  the 
same  way  as  the  air-sea  fluxes  themselves.  The  advantage  of  this  approach  is  it  allows  a  feedback 
from  the  ocean  SST  in  the  heat  flux  calculation  while  at  the  same  time  allowing  for  a  reasonably 
accurate  air-sea  temperature  difference  to  occur.  This  method  will  tend  to  force  the  model  SST  to 
within  1  to  3°C  of  the  prescribed  air  temperature  (Rochford  et  al.  2000),  which  in  turn,  will  reflect 
the  SST  that  was  used  for  the  atmospheric  model  prediction. 

Another  alternative  is  to  use  two-way  coupling  between  the  ocean  model  and  an  atmospheric 
model  so  that  proper  air-sea  feedbacks  between  the  two  models  will  occur  automatically.  This,  of 
course,  is  an  objective  of  the  current  work  with  NCOM.  Such  two-way  coupling  between  atmospheric 
and  ocean  models,  however,  has  a  number  of  issues  of  its  own  (Codron  et  al.  2000). 

7.2  Open  BC 

For  the  surface  elevation  and  normal  barotropic  transport,  the  Flather  radiation  condition  is 
the  most  flexible  BC  to  use.  This  BC  allows  specification  of  the  general  circulation  and  tidal  forcing 
at  the  open  boundary  and  allows  disturbances  within  the  domain  to  radiate  out.  The  clamped  BC 
for  surface  elevation  is  generally  only  useful  for  tidal  forcing  because  it  reflects  outward  propagating 
disturbances  back  into  the  interior,  which  can  disrupt  the  solution.  Even  for  tidal  forcing,  use  of 
a  clamped  BC  requires  a  longer  time  for  the  tidal  solution  to  settle  down  to  a  steady  tidal  cycle 
than  when  using  the  Flather  BC.  This  is  because  of  the  time  required  to  dissipate  transient  tidal 
signals  generated  during  the  spinup  and  reflected  at  the  boundary  by  the  clamped  BC. 

The  advantage  of  the  clamped  BC  for  forcing  tides  is  that  the  normal  velocity  at  the  boundary 
is  not  needed,  i.e. ,  the  BC  is  on  the  elevation  and  the  normal  velocity  is  calculated  by  the  model. 
Because  of  this,  the  clamped  BC  can  be  used  to  calculate  tidal  velocities  at  the  boundary,  which 
can  then  be  used  in  a  more  general  simulation  with  the  Flather  BC.  This  procedure  is  useful  since 
tidal  velocity  data  from  an  external  source  tends  to  be  less  reliable  than  tidal  elevation  data. 

For  the  tangential  barotropic  velocity,  the  simple  zero  gradient  BC  usually  works  fairly  well. 
The  Orlanski  radiation  BC  would  be  a  second  choice.  The  advective  BC  tends  to  respond  too 
slowly  to  changes  at  the  boundary  since  the  advective  time  scale  is  much  slower  than  the  speed  of 
barotropic  signals. 

For  the  normal  baroclinic  velocity,  the  model  calculated  value  works  well,  and  is  generally 
improved  by  the  addition  of  the  upwind  advection  calculation. 
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For  the  tangential  baroclinic  velocity,  as  for  the  tangential  barotropic  velocity,  the  zero  gradient 
BC  seems  to  work  well  and  the  Orlanski  radiation  BC  would  be  a  second  choice. 

For  scalar  fields,  the  Orlanski  radiation  condition  seems  to  work  best. 

8.  SUMMARY 

We  have  described  here  how  BCs  have  been  implemented  in  the  NCOM  for  cyclic  boundaries, 
the  air-sea  interface,  and  along  open  boundaries.  We  have  included  some  explanation  of  the  surface 
BCs  for  the  turbulence  submodels,  and  have  discussed  the  pros  and  cons  of  the  various  BCs  available 
for  the  air-sea  and  open  boundaries.  The  variety  of  formulations  that  can  be  applied  for  the  latter 
gives  NCOM  users  a  great  deal  of  flexibility  in  tailoring  this  OGCM  to  their  available  forcing  func¬ 
tions  and  particular  ocean  applications.  These  can  range  from  one-dimensional  modeling  studies 
with  cyclic  BCs  and  surface  forcing  of  high  temporal  resolution  through  to  basin-scale  and  regional 
coastal  models  forced  by  operational  meteorology  products  from  weather  forecast  centers  and  open 
BCs  from  global  ocean  models.  We  hope  that  this  report  will  be  valuable  to  meteorologists  and 
oceanographers  when  configuring  the  NCOM  for  their  particular  ocean  applications. 
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FORTRAN  SYMBOLS 

The  Fortran  variable  names  used  in  the  NCOM  that  are  cited  within  this  report  are  listed 
below  for  ease  of  reference. 

Al.  Physical  Options 

The  mixing  scheme  used  in  the  NCOM  simulation  determines  whether  open  BCs  will  also  be 
applied  to  the  turbulent  kinetic  energy  (TKE)  and  its  product  with  the  turbulence  length  scale. 

INDZK  type  of  vertical  mixing  to  be  used: 

=1  constant  values  (zkmmin,  zkhmin), 

=2  Mellor-Yamada  level  2  scheme, 

=3  Mellor-Yamada  level  2.5  scheme, 

=4  Mellor-Yamada  level  2.5  scheme  with  option  for  Craig  and  Banner  (1994) 
treatment  of  surface  TKE  flux  and  length  scale. 

A2.  Numerical  Options  and  Parameters 

The  method  that  is  applied  to  solve  for  the  free  surface  mode  affects  how  the  BCs  are  applied 
for  the  elevation  and  barotropic  transport. 

INDBARO  method  of  solution  used  for  free  surface  mode: 

=1  full  explicit,  the  model  is  run  with  a  single  (very  small)  explicit  timestep, 
=2  semi-implicit, 

=3  split-explicit  (not  currently  implemented). 

ASF  Asselin  filter  coefficient  for  temporal  filtering 

EGl,  EG2  temporal  weights  for  surface  elevation  gradient  in  the  barotropic  momentum 
equations  at  the  new  and  central  time  levels. 

A3.  Surface  Forcing  Options  and  Parameters 

Surface  forcing  is  applied  according  to  the  options  set  below.  When  INDSBC  =  0,  all  the 
surface  flux  arrays  are  set  to  zero  and  there  is  no  surface  forcing.  It  overides  all  other  surface 
forcing  flags.  When  INDSBC  =  1,  surface  forcing  is  provided  as  specified  by  the  other  surface 
forcing  flags. 
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INDSBC 

INDATP 

INDTAU 

INDSFT 

INDSFS 

INDSOL 

INDSST 

INDSSS 

RLAXSST 

RLAXSSS 


flag  to  turn  on  or  off  all  surface  forcing: 

=0  no  surface  forcing, 

=1  surface  forcing. 

surface  atmospheric  pressure  forcing: 

=0  no  forcing  is  applied, 

=  1  provided  from  input  data  file, 

—2  provided  from  coupled  atmospheric  model, 
surface  wind  stress  forcing: 

— 0  no  forcing  is  applied, 

— 1  provided  from  input  data  file, 

—2  provided  from  coupled  atmospheric  model, 
temperature  surface  flux: 

=0  no  forcing  is  applied, 

=1  provided  from  input  data  file, 

=2  provided  from  coupled  atmospheric  model, 
salinity  surface  flux: 

=0  no  forcing  is  applied, 

=1  provided  from  input  data  file, 

=2  provided  from  coupled  atmospheric  model, 
solar  radiation  flux: 

=0  no  forcing  is  applied, 

=1  provided  from  input  data  file, 

—2  provided  from  coupled  atmospheric  model. 

sea  surface  temperature  (SST)  relaxed  to  specified  value  via 

a  correction  to  surface  flux  with  rate  defined  by  RLAXSST  (defined  below): 

— 0  SST  relaxation  is  not  used, 

SST  relaxation  is  used. 

sea  surface  salinity  (SSS)  relaxed  to  specified  value  via  correction 
to  surface  flux  with  rate  defined  by  RLAXSSS  (defined  below): 

=0  SSS  relaxation  is  not  used, 

=1  SSS  relaxation  is  used, 

=2  SSS  relaxation  is  used  with  input  file  separate  from  SST  input  file, 
rate  of  relaxation  of  SST  to  specified  value  in  meters/day. 
rate  of  relaxation  of  SSS  to  specified  value  in  meters/day. 


A4.  Lateral  Boundary  Options 

Listed  below  are  the  various  options  for  applying  conditions  for  fields  along  lateral  cyclic  and 
open  boundaries. 

INDCYC  flag  to  denote  periodic  (cyclic)  boundaries: 

—0  no  cyclic  boundaries, 

=4  cyclic  in  x, 

:^5  cyclic  y, 

=6  cyclic  in  x  and  y. 

INDOBC  flag  to  denote  source  of  open  boundary  condition  data: 
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INDOBE 

INDOBR 


INDOBU 


INDOBV 


INDOBVB 


INDTIDE 


=0  no  open  boundaries,  all  open  boundary  points  are  set  to  land, 

=1  use  initial  values  at  boundary  points  for  all  time, 

=2  use  boundary  data  from  input  file, 

=3  use  boundary  data  from  grid  in  which  nested, 
open  boundary  conditions  for  surface  elevation: 

=1  clamped  (elevation  specified), 

=2  Flather  radiation  condition, 

open  boundary  condition  for  scalar  fields  (T  and  S): 

=1  specified  (not  implemented), 

=2  Orlanski  radiation  for  outgoing  and  relaxation  to  specified  value 
for  incoming, 

=3  horizontal  advective  boundary  condition, 

=4  full  3D  advective  boundary  condition. 

open  boundary  condition  for  normal  baroclinic  velocity: 

==1  calculated  by  the  model  ignoring  advection  and  horizontal  mixing, 
=2  Orlanski  radiation  for  outgoing  and  relaxation  to  specified  value 
for  incoming, 

=3  model  calculated  as  in  (1)  plus  upwind  horizontal  advection, 

=4  advective  boundary  condition. 

open  boundary  condition  for  tangential  baroclinic  velocity: 

=1  zero  gradient,  i.e. ,  set  the  same  as  the  interior  value, 

=2  Orlanski  radiation  for  outgoing  and  relaxation  to  specified  value  for 
incoming, 

=3  advective  boundary  condition. 

open  boundary  condition  for  depth-averaged  tangential  velocity: 

=1  zero  gradient,  i.e. ,  set  the  same  as  the  interior  value, 

=2  Orlanski  radiation  condition  for  outgoing  and  relaxation  to  specified 
value  for  incoming, 

~3  advective  boundary  condition. 

include  specified  tidal  forcing  at  open  boundaries: 

=0  no  tides, 

>  0  tides. 


A5.  Lateral  Boundary  Parameters 


These  are  the  time  scales  used  to  relax  to  specified  values  when  the  signal  is  incoming  for  the 
Orlanski  radiation  boundary  condition.  They  are  specified  as  an  e-folding  time  in  model  timesteps. 
Note  that  this  has  the  disadvantage  that  the  results  will  be  dependent  upon  the  time  step  if  the 
values  are  not  changed  when  the  timestep  is  changed. 


RLXOBVB 

RLXOBV 

RLXOBR 


relaxation  time  scale/ for  depth- averaged  velocity, 
relaxation  time  scale/At^  for  baroclinic  velocity, 
relaxation  time  scale/At^  for  scalar  fields. 
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A6.  Lateral  Boundary  Fields 

These  are  the  variable  names  used  to  contain  the  various  fields  along  the  lateral  boundaries. 
Note  that  the  velocities  are  defined  in  terms  of  whether  they  are  normal  or  tangential  to  the 
boundary.  For  example,  the  zonal  baroclinic  component  u  is  normal  at  the  W  and  E  edges  (un) 
and  tangential  at  the  N  and  S  edges  (u^). 

EOB  elevation  at  lateral  boundaries. 

ROB  scalars  at  lateral  boundaries. 

UBOB  barotropic  transport  normal  to  lateral  boundaries. 

For  the  E>-W  boundaries  it  is  U  and  for  the  N-S  boundaries  it  is  V. 

VBOB  barotropic  transport  tangential  to  lateral  boundaries. 

For  the  E-W  boundaries  it  is  V  and  for  the  N-S  boundaries  it  is  U. 

UOB  baroclinic  velocity  normal  to  lateral  boundaries. 

For  the  E-W  boundaries  it  is  u  and  for  the  N-S  boundaries  it  is  u. 

VOB  baroclinic  velocity  tangential  to  lateral  boundaries. 

For  the  E-W  boundaries  it  is  v  and  for  the  N-S  boundaries  it  is  n. 

EOB  elevation  values  at  lateral  boundaries. 
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GRID  INDEXING  FOR  OPEN  BCS 


The  indexing  to  identify  the  locations  of  boundary  grid  points  in  the  NCOM  must  be  robust 
enough  to  handle  open  boundaries,  nesting  of  grids,  domain  decomposition  of  the  grid  into  tiles 
for  parallel  computations,  and  model  domains  with  changing  land  points.  To  accommodate  these 
possibilities,  the  open  boundary  points  are  identified  according  to  a  depth  (or  land-sea)  mask  (H) 
and  an  array  (lEC)  denoting  whether  the  edge  of  a  particular  tile  is  an  exterior  edge  of  the  model 
domain  (lEC  =  1)  or  an  interior  edge  that  abuts  a  neighboring  tile  (lEC  =  0).  Cyclic  boundaries 
(identified  by  INDCYC)  are  not  included  in  the  setting  of  open  boundary  points  as  a  matter  of 
course. 

The  open  boundary  points  are  indexed  in  a  sequential  manner  (NOBPT  =  1,2,...,  NOB)  for 
the  west  (W),  east  (E),  south  (S),  and  north  (N)  boundaries  as  indicated  in  Fig.  Bl.  Indexing  starts 
in  the  lower  left  (SW)  corner  and  proceeds  up  the  left  (W)  edge  to  the  (NW)  corner.  Indexing 
then  continues  from  the  lower  right  (SE)  corner  and  proceeds  up  the  right  (E)  edge  to  the  (NE) 
corner.  Next  it  is  along  the  bottom  (S)  edge  from  the  lower  left  (SW)  corner  to  the  lower  right 
(SE)  corner,  and  finally  along  the  top  (N)  edge  from  the  upper  left  (NW)  corner  to  the  upper  right 
(NE)  corner.  The  index  limits  for  the  open  boundary  points  at  each  of  the  WESN  edges  are  stored 
in  separate  2x4  arrays  for  the  elevation  (NEOB),  the  normal  velocity  (NUOB),  and  the  tangent 
velocity  (NVOB).  Using  elevation  as  an  example,  the  index  limits  for  each  of  the  edges  are  stored 
as  indicated  in  Table  Bl. 

Table  Bl  —  Index  Limits  for  Elevation  Points 


Edge 

Index 

Start 

Index 

End 

W 

NEOB(l,l) 

NEOB(2,l) 

E 

NEOB(l,2) 

NEOB(2,2) 

S 

NEOB(l,3) 

NEOB  (2,3) 

N 

NEOB(l,4) 

NEOB  (2,4) 

To  illustrate  the  indexing  scheme  we  present  in  Table  B2  the  index  limits  for  the  scalars  and 
velocities  for  the  domain  used  in  the  California  Current  System  (CCS)  model  of  the  Coupled 
Biophysical-Dynamics  Across  the  Littoral  Transition  (CoBALT)  project.  The  latter  domain  has 
a  completely  open  western  boundary,  a  completely  closed  eastern  boundary,  and  partially  open 
southern  and  northern  boundaries.  On  the  163  x  229  domain  grid,  the  number  of  open  scalar 
grid  points  is  n  =  229,  0,  162,  and  79,  for  the  WESN  edges,  respectively.  For  u  and  v  along 
the  W  boundary,  the  numbers  of  points  are  227  and  228,  respectively,  because  the  u(l,l),  u(l,m). 
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and  v(l,l)  corner  points  are  not  used  (c.f.  Table  B3).  The  length  of  the  edge  is  defined  as  L  = 
NEOB(2,  *)  -  NE0B(1,  *)  +  1.  Note  that  the  model  program  assigns  the  E  boundary  an  end  index 
value  that  is  less  than  the  starting  index  to  suppress  loop  iteration,  thereby  simplifying  coding. 
For  the  scalar  boundary  points,  the  corner  points  are  always  included  with  the  W  and  E  boundary 
segments  and  not  with  the  N  and  S  boundary  segments. 


Table  B2  —  Index  Limits  for  CCS  Lateral  Boundaries 


NEOB 
Edge  Start  End 

L 

NUOB 

Start  End 

L 

NVOB 

Start  End 

L 

W 

1 

229 

229 

2 

228 

227 

1 

228 

228 

E 

0 

-1 

0 

0 

-1 

0 

0 

-1 

0 

S 

230 

390 

161 

230 

390 

161 

229 

389 

161 

N 

391 

468 

78 

391 

468 

78 

390 

467 

78 

The  index  value  NOBPT  contained  within  the  ranges  NEOB(l,  *)  <  NOBPT  <  NEOB(2,  *) 
provides  the  pointer  to  the  grid  location  (I,  J)  of  the  elevation  boundary  point  via  I  =  lOB(NOBPT) 
and  J  =  JOB(NOBPT),  with  the  same  being  true  for  the  normal  velocity  points  using  NUOB.  The 
indexing  logic  is  the  same  for  NVOB  except  the  grid  locations  are  given  via  I  =  IVOB(NOBPT) 
and  J  =  JVOB(NOBPT).  For  each  elevation  boundary  point,  the  associated  adjacent  interior  grid 
point  (lOBI,  JOBI)  is  identified  as  well  as  its  direction  relative  to  the  boundary  point.  The  grid 
locations  for  the  boundary  points  of  the  tangent  velocity  are  stored  in  separate  arrays  (IVOB, 
JVOB)  because  of  the  separate  indexing  used  for  the  tangent  velocities.  As  an  example,  we  provide 
in  Table  B3  the  (I,  J)  grid  locations  for  the  CCS  lateral  boundary  limits  given  in  Table  B2. 


Table  B3  —  Grid  Locations  (/,  J)  of  CCS  Lateral 
Boundary  Limits 


NEOB 

NUOB 

NVOB 

Edge 

Start 

End 

Start 

End 

Start 

End 

w 

E 

S 

N 

(1,1) 

(1,229) 

(1,2) 

(1,228) 

(1,2) 

(1,229) 

(2,1) 

(2,229) 

(162,1) 

(79,229) 

(2,1) 

(2,229) 

(162,1) 

(79,229) 

(2^1) 

(2,229) 

(162,1) 

(79,229) 

The  indexing  scheme  used  for  the  open  boundary  data  in  the  NCOM  may  or  may  not  seem 
particularly  intuitive.  It  was  chosen  to  save  time  and  storage  for  cases  where  open  boundary  points 
occur  on  only  a  small  part  of  the  external  model  boundary.  Subroutine  OBCPTS  in  the  NCOM 
sets  up  the  locations  and  indexes  of  the  boundary  points  when  a  simulation  is  started  (based  on 
the  model  bathymetry)  and  can  also  be  used  to  define  the  boundary  points  when  setting  up  the 
input  boundary  data  for  a  simulation.  Using  OBCPTS  to  define  the  open  boundary  points  and 
performing  interpolations  in  terms  of  the  location  (latitude,  longitude,  and  depth)  of  the  boundary 
points  can  reduce  the  work  involved  in  setting  up  the  open  boundary  data. 
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For  parallel  processing,  the  NCOM  requires  that  the  overall  horizontal  grid  dimensions  be  an 
integral  multiple  of  the  number  of  tiles  used  in  the  respective  coordinate  direction.  Hence,  in  order 
to  increase  the  flexibility  in  setting  the  model  domain  size,  the  open  boundaries  can  be  inset  from 
the  edge  of  the  grid.  Input  parameter  IBO  specifies  the  number  of  grid  points  that  the  boundary  is 
inset  from  the  edge  of  the  grid  along  each  side  (i.e. ,  the  W,  E,  S,  and  N  sides)  of  the  domain.  If  the 
boundary  is  on  the  edge  of  the  grid,  IBO  is  set  equal  to  zero.  For  the  grid  illustrated  in  Fig.  Bl, 
the  open  boundaries  are  at  the  edge  of  the  grid,  i.e. ,  at  nl  =  1  and  n2  =  n,  ml  =  1,  and  m2  =  m. 
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Pig  31  —  Indexing  used  for  the  NCOM  open  boundaries.  Here  e  denotes  the  elevation,  u 
and  V  the  velocity  components,  B  the  open  boundary,  and  NU  the  grid  points  not  used  in 
the  NCOM. 


Appendix  C 

FINITE  DIFFERENCE  ADVECTIVE  BC  EQUATIONS 


To  make  it  easier  for  the  user  to  understand  the  coding  of  the  advective  BCs,  we  present 
here  the  advective  BC  equations  for  each  of  the  WESN  boundaries  in  finite  difference  form  for  the 
specific  case  of  the  tangential  barotropic  transport.  All  equations  originate  from 


dt  dXn  ’ 


(Cl) 


where  Xn  is  the  coordinate  normal  to  the  boundary  having  increasing  value  in  the  inflow  direction, 
and  inflow  (outflow)  is  the  direction  into  (out  of)  the  model  domain.  The  normal  derivative  of  the 
tangential  component  is  given  by 


dVt  _  j  {Vt-  Vob)/^Xn  >  0  (inflow) 
dXn  ~  \  {Vi-  Vt)/Axn  Un  <  0  (outflow). 


(C2) 


This  differential  equation  has  the  finite  difference  form 

Vt{B,N  +  l)  ^  Vt{B,N)  -  ^[{c+\c\)[VtiB,N)  -  Vob{B,N)]  + 

"  ic~\cm{B,N)-Vt{B,N)]}, 

VniB,N) 
max(e,  D(B,  A))  ’ 


(C3) 


(C4) 


where  Ati  is  the  model  time  step,  Aa:„  is  the  grid  spacing  in  the  normal  direction,  and  e  is  a  small 
constant  to  avoid  division  by  zero.  The  definition  of  Vn  in  Eq.  (C2)  and  Vi  in  Eqs.  (C2)  and  (C3) 
differs  for  each  of  the  WESN  boundaries. 


Cl.  Western  Boundary 

Figure  Cl  shows  the  relevant  variables  for  the  western  boundary.  For  this  boundary,  the 
variables  arc 


i  [U(ivob  +  1,  jvob  -  1)  +  U(ivob  +  1,  jvob)]  (C5) 

Vt  =  V(ivob,jvob)  (C6) 

Vi  =  V(ivob  +  1,  jvob)  (C7) 

Axn  =  Ax.  (C8) 


For  this  case  the  normal  coordinate  is  =  x  and  yields  for  the  normal  derivative 
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dVi  _  j  {Vt-  Vob)/^x  Ki  >  0  (inflow) 
dxn  ~\  {Vi-  Vt)/Ax  Vn<0  (outflow). 


(C9) 


Note  the  open  boundary  lies  along  the  normal  velocity  points,  and  that  the  tangential  velocities  in 
Fig.  Cl  therefore  actually  lie  between  the  boundary  value  (B)  and  the  interior  value  (In). 


Out 

B 

In 

1 

— >  inflow 

X 

X 

X 

1 

outflow 

146 

Vi 

Fig.  Cl  —  Tangential  barotropic  transports  at  the  western  boundary  of  the  NCOM 
domain.  The  boundary  is  denoted  by  B.  The  region  to  the  left  (Out)  is  outside  the 
NCOM  domain  while  the  region  to  the  right  (In)  is  inside.  The  direction  of  inflow/outflow 
is  indicated  by  the  arrows. 

C2.  Eastern  Boundary 


Figure  C2  shows  the  relevant  variables  for  the  eastern  boundary, 
variables  are 

Vn  =  [U(ivob,  jvob  -  1)  +  U(ivob,  jvob)] 

Vt  =  V(ivob,  jvob) 

Y-  =  V(ivob  —  1  Jvob) 

Axn  =  Ax, 


For  this  boundary,  the 

(CIO) 

(Cll) 

(C12) 

(C13) 


Here  the  normal  coordinate  is  Xn  =  —x  and  yields  the  same  expression  for  dVtjdxn  as  in  the  case 
of  the  western  boundary. 
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Fig.  C2  —  Same  as  in  Fig.  Cl,  but  for  the  eastern  boundary 

C3.  Southern  Boundary 

Figure  C3  shows  the  relevant  variables  for  the  southern  boundary.  For  this  boundary,  the 


variables  arc 

Vn  ^  [V(ivob  —  1  Jvob  +  1)  +  V(ivob,  jvob  +  1)]  (C14) 

Vt  —  U(ivob,jvob)  (CIS) 

Vi  =  U(ivob,  jvob  +  1)  (C16) 

Ax^,  Ay.  (C17) 
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In  X  Vi  t  inflow 

-  X  —  Ff 

Out  X  Vob  i  outflow 

Fig.  C3  —  Same  as  in  Fig.  Cl,  but  for  the  southern  boundary 

C4.  Northern  Boundary 

Figure  C4  shows  the  relevant  variables  for  the  northern  boundary.  For  this  boundary,  the 
variables  are 


Vn  =  -^  [V(ivob  -  1,  jvob)  +  V(ivob,  jvob)] 

(C18) 

Vt  =  U(ivob,  jvob) 

(C19) 

Vi  =  U(ivob,  jvob  —  1) 

(C20) 

Axn  =  Ay. 

(C21) 

Out  X  Vob  i  inflow 

—  X  —  Ft 

In  X  Fj  t  outflow 

Fig.  C4  —  Same  as  in  Fig.  Cl,  but  for  the  northern  boundary 


Appendix  D 
GLOSSARY 


The  acronyms  appearing  throughout  this  report  are  listed  below  for  ease  of  reference. 


pseudo  ID 
pseudo  2D 
BC 
CCS 

COAMPS 

CoBALT 


Fortran 

IR 

NCOM 

OGCM 

PWC 

POM 

sss 

SST 

TKE 


one-dimensional  in  the  vertical 
two-dimensional  vertical  plane 
boundary  condition 
California  Current  System 

Coupled  Ocean- Atmosphere  Mesoscale  Prediction  System 
Coupled  Biophysical-Dynamics  Across  the  Littoral  Transi¬ 
tion 

Formula  translator 
surface  solar  IRradiance 
Navy  Coastal  Ocean  Model 
Ocean  General  Girculation  Model 
Pacific  West  Coast 
Princeton  Ocean  Model 
Sea  Surface  Salinity 
Sea  Surface  Temperature 
Turbulent  Kinetic  Energy 
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